/**********************************************************
 *
 *  Copyright (c) 2003  SeikoEpson Inc
 *  All Rights Reserved
 *
 *  File name : f_sinh.c
 *  Function  :
 *        f_sinh function returns hypublic sin's value.
 *        This file is copied from math.lib of 
 *                  CC33 tool(CC33v40.exe + CCPA45.exe, 
 *                  math of ansilib33v25 )
 *
 *  original file's Revision  :
 *      2000/02/10    first release                 M.Igarashi
 *  Revision  :
 *      2003/04/07    IrumaSoft M.Takeishi   1.st design
 *
 **********************************************************/

#include <f_math.h>
#include <f_smcvals.h>

#define	K0   0.336083852466226007246e0f
#define	K1   0.3181332190110303793958e2f
#define	K2   0.9178748106756360272984e3f
#define	K3   0.6282813250833546179346e4f
#define	S0  -0.1292607311299731203424e3f
#define	S1   0.6282813250833546241924e4f


// ALGORITHM
// 1.|x| > 0.5
//     sinh x = ( e^x - e^-x ) / 2
// 2.|x| <= 0.5 ( avoid the cancellation )
//     sinh x = x + x^3 / 3! + x^5 / 5!
//

//  the architecture of float floating point
//
//   0 30         23 22                0  32 bit
//   -----------------------------------
//  | |   exponent  |  fraction         |
//   -----------------------------------
//
//  |               |                   |
//  |    8 bits     |    23 bits        |
//
//         bit    31         sign bit        (  1 bit  )
//             30 - 23      exponent part    (  8 bits )
//             22 -  0      fraction part    ( 23 bits )
//
//


float f_sinh(float sfX){

	long lX;
	unsigned long ulx;
	float sfTemp,sfTemp2,sfX2,sfX3,sfRet;
	
	F_GETW(lX,sfX);				// get raw data
	F_CLR_SIGN( sfX2, sfX );	// clear sign
	
	ulx = lX&0x7fffffff;		// mask sign

	
	// |x| > 0.5
	
	if( ulx > 0x42140000 ){	// check qarea   x > 5
	
		/* Is NaN ? */
		if ( ( ulx & 0x7fc00000 ) == 0x7fc00000 ) {
			errno = EDOM;
			return f_NAN._F;
		}
		sfTemp = f_exp( sfX2 + F_LS );
		if ( (lX&0x80000000) == 0x0 ){			// positive
			sfRet = sfTemp;
			return sfRet;
		}else{									// negative
			sfRet = -sfTemp;
			return sfRet;
		}
	}
	
	if( ulx > 0x3f000000  ){		//check area         ulx > 0.5
		if ( (lX&0x80000000) == 0x0 ){			// positive
			sfRet=f_exp(sfX2+F_LS)-f_exp(-sfX2+F_LS);
			return sfRet;
		}else{									// negative
			sfRet=-(f_exp(sfX2+F_LS)-f_exp(-sfX2+F_LS) );
			return sfRet;
		}
	}
	
	// |dfX| <= 0.5
	
	sfX3=sfX2*sfX2;	
	sfTemp=K0*sfX3+K1;
	sfTemp=sfTemp*sfX3+K2;
	sfTemp=sfTemp*sfX3+K3;
	sfTemp*=sfX2;
	sfTemp2=sfX3+S0;
	sfTemp2=sfTemp2*sfX3+S1;
	sfTemp/=sfTemp2;

	if ((lX&0x80000000)==0x0){			// positive
		//  performance tunning
		//sfRet = sfTemp;
		//return sfRet;
		return sfTemp;
	}else{									// negative
		//  performance tunning
		//sfRet = -sfTemp;
		//return sfRet;
		return -sfTemp;
	}

}
